Rate-dependent slip boundary conditions for simple fluids 



O 

o 

(N 
D 

O 

(N 



O 



o 
o 



> 

O 

o 



I 

O 

o 



X 



Nikolai V. Priezjev 

Department of Mechanical Engineering, Michigan State University, East Lansing, MI 48824 

(Dated: February 6, 2008) 

The dynamic behavior of the slip length in a fluid flow confined between atomically smooth 
surfaces is investigated using molecular dynamics simulations. At weak wall-fluid interactions, the 
slip length increases nonlinearly with the shear rate provided that the liquid/solid interface forms 
incommensurable structures. A gradual transition to the linear rate-dependence is observed upon 
increasing the wall-fluid interaction. We found that the slip length can be well described by a 
function of a single variable that in turn depends on the in-plane structure factor, contact density 
and temperature of the first fiuid layer near the solid wall. Extensive simulations show that this 
formula is valid in a wide range of shear rates and wall-fluid interactions. 

PACS numbers: 83.50.Rp, 47.61.-k, 83.10.Rs, 47.85.1b 



I. INTRODUCTION 

The interest in modeling of fluid flow in confined ge- 
ometries has recently been revived because of the need for 
optimal design of micro- and nanofluidic devices [l| . For 
systems with large surface to volume ratio, the fluid flow 
can be significantly affected by slip at the liquid/solid 
interface. The existence of slip and its degree strongly 
depend on structural and dynamical properties of the 
interface. The most commonly used Navier model for 
the partial slip boundary conditions states that the liq- 
uid slip velocity is proportional to the rate of shear 
normal to the surface. The proportionality coefficient, 
so-called slip length, is defined as an extrapolated dis- 
tance from the wall where the fiuid tangential veloc- 
ity component vanishes. Alternatively, a ratio of fiuid 
viscosity to the slip length determines a friction coeffi- 
cient at the liquid/solid interface, which relates the in- 
terfacial shear stress and fiuid slip velocity. The slip is 
augmented for specially designed superhydropobic sur- 
faces 0, 0] , high polymer weights 0,0 hydrophobic sur- 
faces with t rapp ed nanobubbles [a, 0, Q and high shear 
rates 0, [U [13] . On the other hand surface rough- 
ness Ejl and hydrophilic surfaces [l^ usually lead 
to a reduction of slip. In spite of the long standing inter- 
est in slip behavior, it is not yet clear how the slip length 
depends on the local shear rate and on microscopic pa- 
rameters of the interface. 

In the last two decades, molecular dynamics (MD) 
simulations were used to examine flow boundary condi- 
tions for simp le fluids confined between atomically flat 
walls [il [H H iH, m, m. The boundary con- 
ditions are very sensitive to the wetting properties and 
molecular roughness of the surface, as well as to the liq- 
uid structure near the wall. In general, the slip is en- 
hanced for weak wall-fluid interactions and incommensu- 
rable periodic structures of the surface potential and the 
first fiuid layer. The slip length was found to correlate 
with the degree of the surface induced order in the ad- 
jacent fiuid layer and wall-fluid interaction energy 
Recently, Barrat and Bocquet 0, |2l| have performed a 
detailed analysis based on the Green-Kubo relation for 



the friction coefficient at the interface to derive a scal- 
ing relation for the slip length dependence on density, 
collective diffusion coefficient and structure factor of the 
first fluid layer near the wall at equilibrium. Dynamics of 
the flrst layer of liquid molecules near the wall is closely 
related to the friction of a monolayer of adsorbed parti- 
cles sliding along a solid substrate 0, [2^, . Molecular 
scale corrugations reduce the effective slip length in cases 
of periodic wall roughness [19, 27], chemically patterned 
surfaces [2^, and atomic roughness due to the variable 
size of the wall atoms [1^ . 

While most of the studies have investigated how a vari- 
ation of surface energy and roughness affects boundary 
conditions, the dynamic behavior of the slip length with 
increasing shear rate has received much less attention. 
Difficulties in analysis of the effective slip arise from a 
combination of different factors, such as surface rough- 
ness, wettability and rate-dependency, which produce 
non-equal or even opposite effects on the fiuid fiow near 
the boundary. For example, in non- wetting systems, a re- 
duction of the slip length due to surface roughness might 
be compensated by rate-dependent effects 13]. Thus, the 
understanding of the dynamic behavior of the slip length 
is important for both the interpretation of experimental 
results for flows past rough surfaces 16 , [sT 
fiuid flows in microfluidic channels [31 1. 



and modeling 



Molecular dynamics simulations [9] of simple fluids un- 
dergoing planar shear flow past atomically smooth sur- 
faces have shown that the slip length increases nonlin- 
early with the shear rate, and the usual Navier slip con- 
dition is only valid in a limit of low shear rates. A later 
study [3^ demonstrated that nonlinear boundary condi- 
tions could also describe the flow of complex fluids which 
consist of short polymer chains. A transition from nega- 
tive to positive slip with varying shear rate in Poiseuille 
flows of simple fluids was observed for atomically rough 
hydrophilic surfaces, while at smaller wall-fluid interac- 
tions, the slip length increased approximately linearly 
with shear rate [33]. Experimental studies have also re- 
ported rate-dependent slip for Newtonian liquids in pres- 
sure driven flows in hydrophobic microchannels pl^ and 
thin film drainage in the surface force apparatus p^ . 
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Currently, there is no consensus regarding the functional 
form of the rate-dependent slip length and the existence 
of a shear rate threshold. As a consequence, this pre- 
vents the analysis of more complex systems involving 
combined effects of surface roughness, wettability and 
rate-dependency. 

The focus of this paper is to explore the influence of 
the wall-fluid interaction energy and shear rate on slip 
flow of simple fluids driven by a constant force. We will 
show that for strong wall-fluid interactions the slip length 
increases linearly with the shear rate provided the liq- 
uid/solid interface forms incommensurable structures. A 
gradual transition in rate-dependence of the slip length, 
from linear to highly nonlinear, is observed upon reduc- 
ing the strength of wall-fluid interactions. A detailed 
analysis of the fluid structure near the solid wall shows 
that in a wide range of shear rates and wall-fluid inter- 
actions the slip length can be expressed as a function of 
a single variable that depends on the in-planc structure 
factor, contact density and temperature of the adjacent 
fluid layer. 

This paper is organized as follows. In the next sec- 
tion, we describe details of molecular dynamics simula- 
tions. Predictions from the continuum hydrodynamics 
are briefly summarized in Section IIIII Simulation results 
for the fluid structure and the slip length are presented 
in Section IIVI The summary and conclusions are given 
in the last section. 



II. SIMULATION MODEL 

The computational domain consists of a monoatomic 
fluid confined between two atomistic walls. The fluid 
molecules interact through the pairwise Lennard- Jones 
(LJ) potential 



tion for a fluid molecule of mass m are given by 



(1) 



where e and a represent the energy and length scales of 
the fluid phase. For computational efficiency the cutoff 
distance is set to Tc = 2.5 ct. The LJ wall- fluid interaction 
energy £„f and the length scale (Jwf are measured in units 
of e and a, respectively. In all our simulations, wall atoms 
do not interact with each other and their diameter cTw is 
equal to a. A constant volume accessible to N= 3456 
molecules corresponds to the fluid density p = 0.81a~^. 

The planar Poiseuille flow was generated by a con- 
stant external force in the x direction, which was added 
to the equation of motion for each fluid molecule. The 
heat exchange between the fluid and an external reser- 
voir was regulated by a Langevin thermostat with a ran- 
dom force and a damping term with friction coefficient 
r = 1.0 r^^. This value of the friction coefficient is small 
enough not to influence signiflcantly dynamics of fluid 
molecules (s^ . Issj . The damping term was only applied 
to t he y coordinate to avoid a bias in the flow direc- 
tion 17|. All three components of the equations of mo- 
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where /, is a randomly distributed force with zero mean 
and variance, {fi(0)fj{t)) = 2mkBTTS{t)dij, determined 
from the fluctuation-dissipation relation. The tempera- 
ture of the Langevin thermostat is set to T =1.1 e/ks, 
where fcs is the Boltzmann constant. The equations of 
motion are integrated using the flfth order gear-predictor 
algorithm [s^ with a time step At = 0.002 t, where 
T = ^ma^ je is the characteristic LJ time. 

The upper and lower walls of the cell each consisted 
of 648 atoms distributed between two (111) planes of 
the face-centered cubic (fee) lattice. A flxed wall density 
Pu, = 2.73 corresponds to a nearest-neighbor distance 
d = 0.8cr between equilibrium positions of wall atoms in 
the xy plane. The distance between planes containing 
wall atoms in a contact with the fluid was set to a con- 
stant value of 24.58(7. The dimensions of the cell in the 
xy plane were fixed to 25.03 cr x 7.22 cr. Periodic bound- 
ary conditions were applied along the x and y directions. 

The steady state Poiseuille fiow was induced by a con- 
stant force in the x direction, while both the lower and 
upper walls remained stationary. Initially, fluid molecules 
were uniformly distributed on the centers of the fee lat- 
tice. After an equilibration period of 100 r, an external 
force fx was gradually increased from zero to its flnal 
value corresponding to a steady state flow during next 
10'^ T. After an additional equilibration of about 10^ r, 
fluid velocity proflles were averaged within slices of the 
computational domain of thickness Az = 0.2 a for a time 
interval up to 2 ■ lO^r. Fluid density proflles were com- 
puted within slices of thickness Az = 0.01 a for a time 
period lO^r. 

Fluid density and interaction parameters used in 
this study correspond to the fluid viscosity [i = (2.2 ± 
0.2) ETCT"'^, which was found to be shear rate indepen- 
dent in a range 7 < 0.16 [1, H^l- The upper estimate 
of the Reynolds number (based on the maximum differ- 
ence in fluid velocities at the center and near walls, fluid 
viscosity, and channel width) is Re«10, indicating lami- 
nar flow conditions. 



III. HYDRODYNAMIC PREDICTIONS 

For the planar Poiseuille flow under an externally ap- 
plied force fx in the x direction, the solution of the Navier- 
Stokes equation is described by a parabolic velocity pro- 
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FIG. 1: (Color online) A snapshot of fluid molecules (o) and 
wall atoms (x) coordinates projected on the xz plane for 
fx = and £wf/e = l-l- Particles positions are only shown in 
a range —3 ^ x/a ^ 3 (left). Averaged density profiles for 
£wf/£ = 0.3 (o), e„f/e = l.l (o) and fx = 0.001 e/cr (right). 
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FIG. 2: (Color online) Instantaneous x and y coordinates of 
the fluid molecules (o) in a contact with the lower wall atoms 
( X ) after an equilibration period of IO'^t for fx = 0. Wall-fluid 
interaction energy is fixed to £wf/e = 0.3 (top) and ewf/e=l.l 
(bottom). The fee wall layer is located at 2; = — 12.29 a. 
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where the fluid viscosity fj, is assumed to be shear rate 
independent. The boundary conditions for the fluid 
velocity are prescribed at the confining parallel walls, 
v{—h) = v{h) = Vs- The shear rate at the liquid/solid in- 
terface relates the fluid slip velocity and the slip length 
as follows 
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A quantity of interest for experimental measurements, 
the flow rate, can be evaluated by integrating the fluid 
velocity profile, Eq. ([5]), across the channel width 
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where the second term represents a correction to the flow 
rate due to slip boundary conditions. A relative increase 
in the flow rate due to slip can also be expressed in terms 
of the slip length and the distance between confining 
walls 
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Fluid velocity profiles obtained from MD simulations will 
be compared with the hydrodynamic predictions, Eq. ([5]), 
in the next section. Parameters of the liquid/solid inter- 
face correspond to a flow regime, where the slip length is 
comparable with the channel width, 2/i; and, therefore, 
the flow rate strongly depends on the boundary condi- 
tions. 



IV. RESULTS 



Fluid structure near the walls 



Dynamical and structural properties of a fluid can 
be significantly affected by the presence of a solid sub- 
strate [S^l • A flat solid wall constrains the motion of fluid 
molecules in a normal direction and induces oscillations 
in the fluid density profile. Typically, these density oscil- 
lations gradually decay within a few molecular diameters 
away from the wall. The spatial distribution of molecules 
near the wall becomes non-isotropic and consists of fluid 
layers with thickness of about a molecular diameter. Al- 
though a large amplitude of density oscillations near the 
wall is a signature of a high surface attraction energy, 
the enhanced layering does not necessarily correlate with 
a reduction of the fluid slippage at the interface. For 
example, a highly attractive surface potential free of lat- 
eral corrugations would induce substantial fluid layering, 
which can be interpreted as an infinite slip. In our simu- 
lations, density oscillations relax to a uniform bulk profile 
within 5—6 a away from the wall, see Fig.[T] As expected, 
the higher surface attraction energy causes a more pro- 
nounced fluid layering. 

In general, the fluid layer closest to the flat wall 
has the largest degree of in-plane order characterized 
by the structure factor 5(k) = 1/Ne \ J2j e^^'^^ |^, where 
Tj — (xj,yj) is a two-dimensional vector of a molecule 
position and Ni is the total number of molecules within 
the adjacent layer. Factors affecting the in-plane or- 
der include a correlation between fluid molecules near 
the wall and the energy landscape of the surface poten- 
tial. The degree of the surface induced order depends 
on a mismatch between the wall lattice constant and the 
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FIG. 3: Structure factor S{k,^,ky) in the first fluid layer for 
£wf/£ = 0.3 (top) and ewf/e = l-l (bottom). A force per fluid 
molecule is fx = 0.001 e/a (a), 0.012 e/cr (b), 0.001 e/a (c), and 
0.025 e/o- (d). A small peak appears at the first reciprocal 
lattice vector Gi = (9.04 cr~\ 0). 



FIG. 4: (Color online) Averaged velocity profiles, {v)r/(T, 
for indicated values of an applied force per fiuid molecule. 
Wall-fiuid interaction energy is set to £wf/e = 0.3 (left) and 
£wf/£ = l-l (right). Solid lines represent parabolic fit to the 
data. Dashed lines denote positions of liquid/solid interfaces. 
Vertical axes coincide with fee lattice planes at z = ± 12.29 a. 



nearest-neighbor distance in the adjacent fluid layer. If 
the wall-fluid interaction e^f is comparable with the fluid- 
fluid energy scale e, then commensurable wall-fluid struc- 
tures would typically result in stick boundary conditions, 
or even fluid epitaxial locking, while incommensurable 
structures would likely produce more slippage [13, Sl| • 

In this study, the liquid and solid phases form an in- 
terface with a mismatch between neighboring distance 
within the flrst fluid layer (of about a) and a lattice con- 
stant in the xy plane, d = 0.8a. Figure[2] shows instan- 
taneous snapshots of molecular positions in a fluid layer 
adjacent to the lower wall for the highest (ewf/e = l-l) 
and the lowest (£wf/e = 0.3) wall-fluid interaction ener- 
gies. In the later case, the averaged structure factor 
within the first layer exhibits typical short-range fluid 
ordering characterized by a circular ridge at |k| « 27r/a 
with an amplitude Si « 2.2, see Fig.|31 For ewf/£ = l-l, 
the surface induces higher short-range order, which is 
enhanced along the crystal axes in the xy plane. The 
height of the largest peak in Fig.[3](a) is 5*1 « 4.1. A 
smaller peak of the in-plane structure factor due to peri- 
odic surface potential corresponds to the first reciprocal 
lattice vector d = (9.04 cr"!, 0), see Fig.|3](a) - (d). The 
amplitude of the peak at Gi decreases at larger values of 
the external force. A correlation between surface induced 
order in the first fiuid layer and the slip length will be 
discussed in the next subsection. 



B. Fluid velocity profiles and slip length 

The magnitude of the external force, which is required 
to reach a parabolic velocity profile described by Eq. ^ , 
depends on the fluid viscosity, density and wall-fluid in- 



teraction parameters [H, [13, [12, [H, [1^ ■ Since the slip 
velocity is not known a priori, the value of the force in 
MD simulations is usually adjusted so that a fluid veloc- 
ity profile could be accurately resolved without an exces- 
sive computational effort due to thermal averaging. One 
of the goals of this study is to systematically explore the 
effect of an applied force on a flow of simple fluids near 
solid boundaries and to determine a variation of the slip 
length as a function of shear rate. In our simulations, 
the channel width, 2h — 23.58 a, is large enough to avoid 
extreme confinement conditions (ssl . [sol . |40| , where devi- 
ations from macroscopic hydrodynamics are expected. 

Examples of averaged velocity profiles for different val- 
ues of the external force fx and fixed wall-fiuid interac- 
tion energies, £wf/£ = 0.3 and e^{/e = l.l, are shown in 
Fig-Hi The data are presented only in half of the chan- 
nel because of the symmetry with respect to ^ = plane, 
v{z) = v(— z). Fluid velocity profiles are well fitted by 
a parabola as expected from the hydrodynamic predic- 
tions, see Eq. ([5]). Weak oscillations within 2 a near the 
walls correspond to a pronounced fluid layering perpen- 
dicular to the surface. In a range of wall-fluid interaction 
energies considered in this study, 0.3 ^Cwf/ 1.1, fluid 
flow undergoes slippage at the solid walls. Fluid velocity 
in the channel and at the interfaces increases with the 
applied force. The shear viscosity, /i= (2.2 ± 0.2) eTa~^, 
which was computed from the Kirkwood relation [4]| . 
remained independent of the applied force. 

A ratio between fiuid slip velocity and the local shear 
rate at the interface defines the slip length, denoted by 
Ls throughout. For the parabolic profiles, the slip length 
was evaluated by linear extrapolation of the slope at the 
interface to zero velocity. In our simulations, the posi- 
tion of the interface in the z direction is defined at a 
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FIG. 5: (Color online) Variation of the slip length, L^/a, as 
a function of an applied force per fluid molecule. Wall-fluid 
interaction energy is set to Ewf/e ~ 0.3 (a), 0.4 (<), 0.5 (o), 0.7 
(o), and 1.1 (>), respectively. Dashed curve is the best fit to 
L.(fx) = L° (l-fx/fc)'°-^ with L° = 19.5a and ic = 0.021 e/a. 
Straight dashed lines represent the best fit to the data. 



FIG. 6: (Color online) Behavior of the slip length as a function 
of the local shear rate at the interface. Values of the wall- fluid 
interaction energy are listed in the inset. The same data as in 
Fig.[5l Dashed curve is the best fit to Eq. ^ with L° = 19.5 a 
and 7c = 0.093 r~^. Solid curves are a guide for the eye. 
Dashed lines show the best linear fit to the data. 



distance 0.5 CTw away from the fee lattice planes, see ver- 
tical dashed lines in Fig.[4l This offset was chosen to 
account for the excluded volume due to wall atoms. The 
slip length was defined as the average of values extracted 
at the top and the bottom walls. 

The dynamic response of the slip length as a function 
of the external force is presented in Fig. [51 A gradual 
transition in the functional dependence of is (fx) is ob- 
served by varying the strength of the wall-fluid interac- 
tion. The slip length increases monotonically with the 
applied force for £wf/e5^0.7. The data can be well fitted 
by a straight line, see Fig.[5l At lower surface energies, 
£wf/e^0.5, the relation between Lg and fx becomes non- 
linear and exhibits a pronounced downward curvature 
for Ewf/e = 0.3. For each curve shown in Fig.O a ratio of 
maximum slip length to its value at small applied forces 
is equal to 1.63 ± 0.13. This factor determines an upper 
bound for the increase in the fiow rate due to slip de- 
pendence on the applied force. Given the fixed channel 
width used in our study, the maximum relative gain in 
the fiow rate due to variation of the slip length with the 
force is 1.59 ± 0.08 for £„f/e = 0.3 and 1.36 ± 0.08 for 
£wf/e = l-l- These results suggest that a significant drag 
reduction for laminar flows can be achieved through the 
increase in the pressure difference across a microfiuidic 
channel. 

In the original paper by Thompson and Troian [9] on 
shear flow of simple fluids, the slip length was found to 
increase nonlinearly with the shear rate. In a range of 
accessible shear rates and weak wall-fluid interactions, 
the MD data were well described by a power law function 

L,(7)=L: (1-7/7,)-0-5, (9) 

where L° and 7c are fitting parameters. In our simula- 



tions, the shear rate is proportional to the external force, 
see Eq. ([5]), and, therefore, the analogous expression for 
the slip length dependence on the applied force should 
be L^(fx) = L° (1 - fx/fc)"°-^. This form was used to 
fit the data for the lowest wall-fluid interaction energy 
Ewf /e = 0.3, see dashed curve in Fig.O The agreement is 
rather good in the range of applied forces fx/ fc^ 0.67. 

For each curve presented in Fig. [51 the external force 
was varied from fx = 0.001 e/cr up to a maximum value, 
which depends on e^f- This value of the force corre- 
sponds to a maximum shear stress the liquid/solid inter- 
face can support. With a further increase of the force, 
the fluid flow acquires a large velocity component in the 
X direction, (v) ^ vt, where v'^ = kBT/m is the thermal 
fluid velocity. In this extreme regime, the dynamics of 
fiuid molecules near walls cannot be resolved accurately 
with the integration time step used in this study. We 
note, however, that test runs with a smaller time step 
At = 0.001 T did not produce noticeable changes in the 
results presented in Fig. [51 The transition to the fiow 
regime characterized by very large slip velocities is not 
the main focus of this paper and, therefore, it was not 
studied in detail. 

The parabolic shape of the fluid velocity profiles im- 
plies that the external force fx is proportional to the in- 
terfacial shear rate, see Eq. ([5]). The functional depen- 
dence of the slip length, therefore, is expected to be sim- 
ilar for fx and the local shear rate. FigurelHl shows the 
same MD data as in Fig. [51 but replotted in the axis Lg 
versus local shear rate, as extracted from parabolic ve- 
locity profiles at the location of interfaces. The range of 
shear rates is below the values reported for laminar flows 
in Ref. [l3|. The slip length increases with shear rate, 
and the growth of Lg is enhanced at lower values of Swf ■ 
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FIG. 7: (Color online) A correlation between the slip length, 
Ls/(J, and the inverse value of the in-plane structure factor, 
1/5(01), evaluated at the first reciprocal lattice vector. Solid 
curves are a guide for the eye. The inset shows the same data 
plotted as a function of [S(Gi) pc 
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FIG. 8: (Color online) Log-log plot of the slip length as a 
function of a combined ratio Tiks/e [^(Gi) pcO"^]~^. Values 
of the wall-fiuid interaction energy are tabulated in the inset. 
The same data as in Fig. [5] The solid line with a slope 1.44 
is plotted for reference. 



The dependence ^5(7) for e^f/e = 0.3 can be well fitted 
by the power law function, Eq. ([5]), with L° — 19.5 a and 
7c = 0.093 T~^, see dashed curve in Fig. [HI Thus, for a 
weak wall-fluid interaction our results are in agreement 
with those reported in a previous study ^ on dynamical 
behavior of the slip length in boundary driven flows. Fur- 
thermore, at higher surface energies, £wf/e^0.7, the slip 
length increases linearly with the interfacial shear rate, 
see Fig.ini In this regime, the simulations results might 
be relevant to a monotonic growth of the slip length with 
the shear rate measured in pressure driven flows in hy- 
drophobic microchannels ^12J. Finally, in contrast to our 
results, a transition from the linear to the power law rate- 
dependence upon increasing the strength of wall-fluid in- 
teractions was reported in Ref. [s^ . The difference in the 
slip behavior might be explained by the lower fluid den- 
sity and higher shear rates examined in that study (ssj . 

Molecular-scale corrugations of a solid wall composed 
of periodically arranged LJ atoms induce in-plane order 
in the adjacent fluid layer. The amount of surface in- 
duced order in the first fluid layer is reflected in the fluid 
in-plane structure factor. A correlation between the slip 
length in a shear rate independent regime and the peak 
value of the fluid structure factor evaluated at the first 
reciprocal lattice vector was established earlier [l3, [2l| . 
In this study, the peak of the structure factor at the 
first reciprocal lattice vector Gi = (9.04 ct^^, 0) is dis- 
placed from the circular ridge at the vector |k| « 27r/fT, 
see Fig. [3] The magnitude of the peak at Gi decreases 
with increasing slip velocity. Figure [7] shows the behav- 
ior of the slip length as a function of the inverse value 
of the structure factor. For the largest wall-fluid interac- 
tion energy eyjf/e = l.l the slip length increases linearly 
with l/S'(Gi). At lower surface energies, the function 



deviates from the linear dependence and increases more 
rapidly for ewf/e = 0.3. These results demonstrate that 
the slip length in a shear rate dependent regime strongly 
correlates with the surface induced order in the fluid layer 
adjacent to the wall. 

In a zero shear limit, the Green-Kubo analysis for the 
friction coefficient at the liquid/solid interface shows that 
for attractive wall-fluid interactions the slip length de- 
pends on the structure factor, contact density, diffusion 
coefficient and temperature of the first fluid layer [2l| . In 
our simulations, the equilibrium fluid density and tem- 
perature profiles near the wall are modified at higher 
shear rates. The contact density, pc^ was determined 
from the maximum of the fluid density profile in the first 
fiuid layer, see Fig.[TJ The contact density increases with 
the strength of wall-fluid interaction and decreases with 
shear rate. In the inset of Fig. [7] the slip length is plotted 
against an inverse product of the structure factor and the 
contact density. Except for the lowest wall-fluid interac- 
tion energy, ewf/£ = 0.3, the functional form of the slip 
length consists of nearly linear interconnected segments 
each characterized by its own value of £wf . 

The data for the slip length for different wall-fiuid in- 
teraction energies and shear rates can be collapsed onto 
a single master curve by taking into account a variation 
in temperature of the first fiuid layer. In a steady state 
fiow induced by a constant force, fiuid temperature was 
computed from the kinetic energy 

N 

kBT^^Y.'^v,-^{yi)]\ (10) 
i— 1 

where is a three-dimensional vector of a molecule po- 
sition and v(rj) is the local average fiow velocity. At 
low shear rates, 7 < 0.02 t~^, the fluid temperature re- 
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mains equal to that imposed by the Langevin thermo- 
stat, T =l.le/kB, and it increases shghtly at higher 7. 
The maximum relative increase of T is about 3.5% for 
each value of Swf- The heatup is larger near the walls 
because of the higher shear rates and the slip velocity, 
which becomes comparable to the thermal fluid velocity, 
see Fig.m Temperature of the first fluid layer, Ti, rises by 
about 10% at the highest shear rates reported in Fig. [HI 
Figure [8] shows the slip length as a function of a com- 
bined ratio of temperature to the product S{Gi)pc eval- 
uated in the first fluid layer. In a wide range of shear 
rates and for 0.3 ^ Ewf/e ^1-1, the slip length is well 
described by a power law function 

Ls ^ (Ti/5(Gi)p,)", (11) 

with a = 1.44±0.10. A solid line in Fig. [8] corresponds to 
the value a — 1.44. This result implies that the condition 
Ti/ S{Gi)pc — const defines a curve on the plane £wf and 
7 characterized by a constant slip length. The functional 
dependence of the slip length, Eq. pTjl . can also be deter- 
mined from equilibrium measurements of {S (Gi) pc)~" 
for different surface energies £wf- This prediction is in 
qualitative agreement with previous MD results [13, [HI , 
which demonstrated that the slip length decreases with 
the contact density and surface induced order in the ad- 
jacent fluid layer. A strong correlation between the slip 
length and microscopic properties of the first fluid layer 
provides a framework for the analysis of systems with 
combined effects of wettability and rate-dependency. 

V. CONCLUSIONS 

In this paper the effect of surface energy and shear rate 
on the slip length in a flow of simple fluids was studied 
by molecular dynamics simulations. Fluid velocity pro- 
files in steady state flow induced by a constant force were 



fitted by a parabola with a shift by the value of the slip 
velocity. The slope of the parabolic fit at the interface 
was used to define the local shear rate. For a weak wall- 
fiuid interaction the slip length increases nonlinearly with 
the shear rate and its dependence can be well fitted by a 
power law function. Increasing the strength of wall-fluid 
interactions leads to the linear rate-dependence of the slip 
length. For a fixed channel width, the flow rate increases 
significantly due to the rate-dependence of the slip length 
for both weak and strong wall-fluid interactions. Simu- 
lation results also indicate a strong correlation between 
the slip length and the surface induced order in the first 
fiuid layer in a shear rate-dependent regime. We showed 
that in a wide range of wall-fluid interaction energies and 
shear rates the slip length is well described by a function 
of a single variable that depends on the in-plane struc- 
ture factor, contact density and temperature of the first 
fiuid layer. 

Future work will show how sensitive these results are 
to the variation of molecular-scale roughness. The sur- 
face induced order in the adjacent fluid layer and the 
slip length might be affected by the presence of substrate 
inhomogeneities. The effect of thermal, random and pe- 
riodically corrugated surfaces on slip behavior in a rate- 
dependent regime should be explored. 
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